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Abstract 

The two-site Bose-Hubbard model is a simple model used to study Josephson 
tunneling between two Bose-Einstein condensates. In this work we give an overview 
of some mathematical aspects of this model. Using a classical analysis, we study the 
equations of motion and the level curves of the Hamiltonian. Then, the quantum 
dynamics of the model is investigated using direct diagonalisation of the Hamilto- 
nian. In both of these analyses, the existence of a threshold coupling between a 
delocalised and a self-trapped phase is evident, in qualitative agreement with ex- 
periments. We end with a discussion of the exact solvability of the model via the 
algebraic Bethe ansatz. 
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1 Introduction 



The phenomenon of Bose-Einstein condensation, while predicted long ago [1,2], is nowa- 
days responsible for many current perspectives on the potential applications of quantum 
behaviour in mesoscopic systems. This point of view has arisen with the experimental 
observation of condensation in systems of ultracold dilute alkali gases, realised by several 
research groups using magnetic traps with various kinds of confining geometries [3,4]. 
These types of experimental apparatus open up the possibility for studying quantum 
effects, such as Josephson tunneling and self-trapping [5,6], in a macroscopic setting. 

From the theoretical point of view, the two-site Bose-Hubbard model (see eq. (fT|) 
below), also known as the canonical Josephson Hamiltonian [7], has been a useful model 
in understanding tunneling phenomena. The simplicity of the model means that it is 
amenable to detailed mathematical analysis, as we will discuss below. However despite 
this apparent simplicity, the Hamiltonian captures the essence of competing linear and 
non-linear interactions, leading to interesting, non-trivial behaviour. 

The Hamiltonian is given by 

H = ^(iV 1 -iV 2 ) 2 -^(iV 1 -iV 2 )-^(6t 1 6 2 + 46 1 )- (1) 

where b\,b\ denote the single-particle creation operators in the two wells and N% = 
b\bi,N2 = b\p2 are the corresponding boson number operators. The total boson num- 
ber Ni + N 2 is conserved and set to the fixed value of N. The coupling k provides the 
strength of the scattering interaction between bosons, Afj, is the external potential and 
Sj is the coupling for the tunneling. The change Sj — > —£j corresponds to the unitary 
transformation b\ — > bi, 6 2 — > —bi-, while A/x — > —A/i corresponds to b% <-> 6 2 . Therefore 
we will restrict our analysis to the case of £j, A/i > 0. For k > 0, following [7] it is useful 
to divide the parameter space into three regimes; viz. Rabi (k/Sj << N^ 1 ), Josephson 
(iV -1 << k/Sj << N) and Fock (N << k/Sj). For these three regimes, there is a cor- 
respondence between (JTJ and the motion of a pendulum [7]. In the Rabi and Josephson 
regimes this motion is semiclassical, in contrast to the Fock regime. For both the Fock 
and Josephson regimes the analogy corresponds to a pendulum with fixed length, while 
in the Rabi regime the length varies. 

In the present work, we give an overview of some of the mathematical aspects of 
(JTJ). We undertake an analysis of the classical and the quantum dynamics of the system, 
and discuss how the system exhibits a threshold coupling, originally identified in [8], 
about which the dynamics abruptly changes. Below this threshold point the dynamics 
is delocalised, while above it the dynamics turns out to be localised (macroscopic self- 
trapping). This can be seen at both the classical and quantum level. The result is in 
qualitative agreement with experimental results [6]. We conclude by giving an outline of 
the algebraic Bethe ansatz solution of (JTJ. 
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2 Classical dynamics 

First we study a classical analogue of the model. Let Nj, 9j, j = 1, 2 be quantum variables 
satisfying the canonical relations 

[9x, 6 2 ] = [N h N 2 ] = 0, [N j} 9 k ] = i5 jk I. 

Using the fact that 

ex.p(i9j)Nj = (Nj + l)exp(i^) 
we make a change of variables from the operators bj, bp j = 1, 2 via 

bj = exp(i9j) y/Nj, &j = \/~Njexp(—i6j) 

such that the Heisenberg canonical commutation relations are preserved. Now define the 
variables 

z = (N 1 - N 2 )/N 
<P = N{4> x - 2 )/2 

where z represents the fractional occupation difference (or the imbalance) and (f) the phase 
difference. In the classical limit where N is large, but still finite, we may equivalently 
consider the Hamiltonian [9,10] 

H(z, 0) = {^z 2 -f3z- v / r^cos(20/iV)) (2) 

where 

kN Afi 

and (z, (p) are canonically conjugate variables. We note the Hamiltonian (j2J) obeys the 
symmetries 

H{z,4>)\^ = -H{z,<t> + m/2)U^ 

H{zMxfi = H{-z,4>)\ K _ p . (3) 
The classical dynamics is given by Hamilton's equations of motion 

* = f -¥(^-^7T^ cosW/JV) 



dH 



-£j (yi - z 2 sin(20/AO) . (4) 



Now we study the fixed points of the Hamiltonian (J2J), determined by the condition 
z = (p = 0. This leads to the following classification: 



and z is a solution of 



\z-f3 = (5) 



which has a unique real solution for A > 0. 
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• = Ntc/2 and z is a solution of 

\z-(3 = —L=. (6) 

This equation has either one, two or three real solutions for A > 0. 

From eq. © we can determine that there are fixed point bifurcations for certain 
choices of the coupling parameters. These bifurcations allow us to divide the coupling 
parameter space in two regions. Setting f(z) = Xz — (3 and g(z) = z(l — z 2 )' 1 ^ 2 , the 
boundary between the regions occurs when f(z) is the tangent line to g(z) at some value 
z . A standard analysis shows this occurs when A = g'(z ) = (1 — Zq)~ 3 / 2 . Requiring 
f{ z o) — g( z o) then yields the following relationship 

A = (1 + \(3\ 2/ Y 2 (7) 
determining the boundary. This is depicted in Fig. 




Figure 1: Coupling parameter space diagram identifying the different types of solutions 
for equation (|5]l. In region I there is just one solution for z, a local maximum. In region 
II there are three solutions for z, two local maxima and a saddle point. The boundary 
separating regions I and II is given by (J7J). 

This leads us to the following classification: 

• < A < 1: For any value of (3 there is just one real solution, for which the Hamil- 
tonian attains a local maximum. 

• A > 1: Here transition couplings ±/3q appear, which can be seen from Fig. ^ For 
(3 G (—Po,Po), the equation has two locally maximal fixed points and one saddle 
point, while for (3 > (3q or (3 < —(3q the equation has just one real solution, a locally 
maximal fixed point. 
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We remark that in the absence of the external potential (A/i = (3 = 0) the transition 
value is given by Ao = 1. Using the symmetry relation Q we can deduce that for the 
attractive case A < 0, Ao = — 1 is the coupling marking a bifurcation between a locally 
minimal fixed point (for A > — 1) and two locally minimal fixed points and a saddle point 
(for A < —1). This is a supercritical pitchfork bifurcation of the classical ground state. 
The results of [14] predict that the ground-state entanglement, as measured by the von 
Neumann entropy, is maximal at this coupling. The numerical results of [15] confirm this. 

Next we look at the dynamical evolution. In that which follows we will consider the 
equations (JH) in the absence of the external field (A/z = or, equivalently, (3 = 0). 
An analysis including the effect of this term can be found in references [11, 12]. We 
integrate to find the time evolution for the imbalance z, using the initial condition 
z(0) = 1, 0(0) = 0. By plotting z against the time, it is evident that there is a threshold 
coupling A c = 2 separating two different behaviours in the classical dynamics, as can be 
seen in Fig. |2J 

(i) For A < 2 the system oscillates between z = — 1 and z = 1. Here the evolution is 
delocalised; 

(ii) For A > 2 the system oscillates between z = and z = 1. Here the evolution is 
localised. 

The threshold occuring at A c = 2 was first observed in [8]. 




10 20 30 
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Figure 2: Time evolution for the imbalance z. The solid line is for A = 1.9, while the 
dashed curve is for A = 2.1. Here we are using N = 100, £ j = 1 and the initial conditions 
z(0) = 1, (f)(0) = 0. The threshold coupling occurs at A c = 2. 
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To help visualise the classical dynamics, it is useful to plot the level curves (con- 
stant energy curves) of the Hamiltonian (J2J) in phase space. Given an initial condition 
(z(0), 0(0)), the system follows a trajectory along the level curve H(z(0), 0(0)). In Fig. |3] 
we plot the level curves for different values of A (A = 1.5 on the left and A = 2.5 on the 
right), where we take 20/iV G [— 7r, 7r]. We can observe clearly two distinct scenarios: 

• A > 2: Here we see that for the orbit with initial condition zq = 1, 0(0) = 0, 
increases monotonically (running phase mode). The evolution of z is bounded in 
the interval [0, 1], leading to localisation (self-trapping). 

• A < 2: Here we see that for the orbit with initial condition z(0) = 1, 0(0) = 0, 
the evolution of is oscillatory and bounded in the interval (—Nn/2, Nix/2). The 
evolution of z is not bounded, leading to delocalisation. 



(a) (b) 



— 71 o 71 — 71 q 71 




-71 71-71 71 

24>/N 2t)l N 



Figure 3: Level curves of the Hamiltonian (J2J) (a) for A = 1.5 (below the threshold point) 
and (b) for A = 2.5 (above the threshold point). We are using N = 100 and Sj = 1. 
Above the threshold coupling running phase modes occur leading to localised evolution 
of z. Below the threshold coupling the evolution of z is delocalised. 

The threshold coupling A c = 2 (or k/Sj = 4/iV, in terms of the original variables) 
separates two distinct dynamical behaviours. This value for the threshold between delo- 
calisation and self-trapping also occurs for the quantum dynamics, as we will show in the 
next section. 

3 Quantum dynamics 

We will investigate the quantum dynamics of the Hamiltonian in the absence of the 
external potential (A/j, = 0) using the exact diagonalisation method. It is well known 
that the time evolution of any state is determined by = U(t)\(po), where U is 
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the temporal evolution operator given by U(t) = ^2 m=0 \m)(m\ exp(— iE m t), \m) is an 
eigenstate with energy E m and |0o) represents the initial state. Using these expressions 
we can compute the expectation value of the relative number of particles 

((JVi - N 2 )(t)} = (tf (t)^ - N 2 \V(t)). (8) 

From Fig. H]it is seen that the qualitative behaviour in each region does not depend on the 
number of particles. We find that in the interval k/Sj G [1/N 2 , 1/N] (close to the Rabi 
regime) the collapse and revival time takes the constant value t cr = An 1 In the interval 
between k/Sj = 1/N and k/Sj = 1 the system undergoes a transition from oscillations 
which vary between positive and negative values of (Ni — N 2 ) (delocalised) to one where 
(Ni — N 2 ) is close to N (self-trapping). 



N = lOO 
N = 400 




0.999 I 1 1 1 1 1 1 1 1 1 

O 0.2 0.4 0.6 O.S 1 

t 

Figure 4: Time evolution of the expectation value for the relative number of particles 
for different ratios of the coupling k/Sj from the top (Rabi regime) to the bottom (Fock 
regime): k/Sj = 1/N 2 , 1/N, 1, N, N 2 for N = 100,400 and the initial state is \N,0). 

Now we focus in more detail the time evolution of the expectation value of the relative 
number of particles in the interval k/Sj G [1/N, 1] In Fig. El we present the case iV = 100: 
we observe the evolution of the dynamics from a collapse and revival sequence for k/Sj < 
4/iV, through the self-trapping transition at k/Sj = 4/iV, and toward small amplitude 
harmonic oscillations in the imbalance of the localised state when k/Sj = 1. It is also 
interesting to observe in the localised phase k/Sj > 4/N the re-emergence of a collapse 
and revival sequence. Further increases in k/Sj lead to a decaying of the collapse and 

1 The ratio k/£j = 1/N 2 means that we are using k = 1 and £j = N 2 and similarly for the other 
cases. 



7 



revival sequence toward harmonic oscillations which occur at k/Bj = 1. A more detailed 
investigation, using another initial conditions can be found in ref. [13]. 




t t 

Figure 5: Time evolution of the expectation value between k/Sj = 1/N and k/Sj = 1. 
On the left, from the top to the bottom k/Ej = 1/N,2/N,3/N,4/N and on the right, 
from the top to the bottom k/Sj = 5/N, 10/N,50/N, 1, where N = 100 and the initial 
state is \N, 0). 



From the above picture it is clear that the threshold coupling k/Sj = 4/N predicted by 
the classical analysis, representing the boundary between a delocalised evolution (k/£j < 
4/JV) and self-trapped evolution (k/Sj > 4/N), also holds for the quantum dynamics. 



4 Exact Bethe ansatz solution 



In this final section we briefly discuss the exact Bethe ansatz solution of (0). More details 
can be found in [16, 17]. We begin with the 5'?7(2)-invariant .R-matrix, depending on the 
spectral parameter u: 

( 1 \ 

b(u) c{u) 
c(u) b{u) 
\ 1 / 

with b(u) = u/{u + 77) and c(u) — rj/(u + 77). Above, 77 is an arbitrary parameter, to be 
chosen later. It is easy to check that R(u) satisfies the Yang-Baxter equation 



R{u) 



(9) 



R 12 {u - v)R 13 (u)R 2 3(v) = R23(v)R 13 (u)R 12 (u - v). 



(10) 
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Here Rjk{u) denotes the matrix acting non-trivially on the j-th and k-th spaces and as 
the identity on the remaining space. 

Next we define the Yang-Baxter algebra T(u), 



C(u) D(u) 
subject to the constraint 

Rab(u - v)T a (u)T h (v) = T b (v)T a (u)R ab (u - v). (12) 
We may choose the following realization for the Yang-Baxter algebra 

7c(T a (u)) = L al (u + u)L a2 {u-u), (13) 
written in terms of the Lax operators [16] 

L,(«)=("+r< ,= 1,2. (14) 

Since L(u) satisfies the relation 

Rab(u ~ v)L ai (u)L bi (v) = L bi (v)L ai (u)Ri 2 {u -v), i = 1, 2 (15) 

it is easy to check that the Yang-Baxter algebra ()12j) is also obeyed. 
Finally, defining the transfer matrix through 

t(u) = 7r(ti a T a (u)) = tt(A(u) + D(u)) (16) 

it follows from (|T2*|) that the transfer matrix commutes for different values of the spectral 
parameters; i.e., the model is integrable. Now it is straightforward to check that the 
Hamiltonian (0) is related with the transfer matrix t ()16|) through 

H — —k (t(u) - ^(t'(0)) 2 - wt'(O) - r]- 2 + u 2 -u 2 

where the following identification has been made for the coupling constants 

— = , = —K71UJ, = K. (17) 

4 2 2 ' 2 v ; 

We can apply the algebraic Bethe ansatz method, using the Fock vacuum as the 
pseudovacuum, to find the Bethe ansatz equations (BAE) 

N 

V 2 (v 2 -u J 2 ) = H Vl ~ V \ r] (18) 
and the energies of the Hamiltonian (see for example [16, 17]) 



N 



,2 7\r2 



E = -k[ V - 2 Y\(1 + — !—)-^^-u V N-u 2 

N \ 

- V - 2 + u 2 + {u 2 - u 2 )\{{l - . (19) 



i=l 



Surprisingly this expression is independent of the spectral parameter u, which can be 
chosen arbitrarily. 

Using the Bethe ansatz solution, it is possible to derive form factors and correlation 
functions. Details are given in [16]. 
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